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Abstract 

During the 2010 rainy season in Yangbajing (4300 m above sea level) in Tibet, China, a long- 
duration count enhancement associated with thunderclouds was detected by a solar neutron tele- 
scope and neutron monitors installed at the Yangbajing Comic Ray Observatory. The event, lasting 
for ~40 min, was observed on July 22, 2010. The solar neutron telescope detected significant 7-ray 
signals with energies >40 MeV in the event. Such a prolonged high-energy event has never been 
observed in association with thunderclouds, clearly suggesting that electron acceleration lasts for 40 
min in thunderclouds. In addition, Monte Carlo simulations showed that > 10- MeV 7 rays largely 
contribute to the neutron monitor signals, while >l-keV neutrons produced via a photonuclear 
reaction contribute relatively less to the signals. This result suggests that enhancements of neu- 
tron monitors during thunderstorms are not necessarily a clear evidence for neutron production, 
as previously thought. 

PACS numbers: 52.38.Ph,82.33.Xj, 92.60.Pw, 93.30.Db 
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I. INTRODUCTION 



Recent observations have shown that thunderclouds are powerful particle accelerators, 
emitting bremsstrahlung 7 rays that extend to 10 MeV or higher IK 7 ]]. Unlike terrestrial 
gamma ray flashes (TGFs) and lightning-related 7 rays, which last for milliseconds or less, 
thundercloud-related 7 rays are characterized by durations of a few tens of seconds to a few 
minutes, or occasionally more than 10 min. These thundercloud-related 7 rays have been 
thought to be produced by relativistic electrons, in accordance with the relativistic runaway 
electron avalanche (RREA) model 10] that involves acceleration and multiplication of 
ambient electrons. However, mainly because of the lack of a large sample of thundercloud- 
related 7 rays, there is still no consensus that all of those 7 rays are really generated by the 
RREA mechanism. It is also unclear whether the charging mechanism of thunderclouds is 
related to the production of thundercloud-related 7 rays. 

Several groups conducting their experiments on high mountains have reported detect- 



ing various particles besides 7 rays, in possible association with thunderstorms 11 
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Among such particles, the production of neutrons in coincidence with natural lightning 
by a thermonuclear reaction 2 H + 2 H — > n (2.45MeV ) + 3 He was closely investigated in the 
1970s— 1980s because such neutrons provide a key not only to elucidate the mechanism of 
lightning but also to know if such neutrons are another source captured by 14 C. Investigating 



such neutrons, Shah et al. 



131 ] and Shyam and Kaushik 14j reported detections of 10 7 — 10 



110 



neutrons per lightning strike. 

Instead of the above fusion mechanism, Bahich and Roussel-Dupre 
tonuclear reaction, 14 N(7, n) 13 N, showing that the fusion mechanism is not feasible under the 



15] proposed a pho- 



usual physica 
of 10.5 MeV 



conditions in lightning. The photonuclear reaction begins at a 7-ray energy 



161 ]. and hence may occur because 7 rays with energies above the threshold 



have been actually observed. Therefore, photonuclear neutrons provide another clue to solve 
non-thermal mechanism in thunderstorms. Actually, Carlson et al. jl3] made a close inves- 
tigation on neutron production in TGFs, predicting that a TGF averagely produces ~10 12 
neutrons corresponding; to a ground-level neutron fluence of (0.03 — l)xl0 4 m -2 . Similarly, 
Babich et al. 18] also predict that neutrons with fluence of 10 3 — 10 7 m~ 2 would arrive at 
ground level when energetic 7 rays are produced under the RREA mechanism. Many neu- 
tron monitors, installed at cosmic-ray stations in the world, could detect such neutron bursts 



in thunderstorms, if those neutrons actually reaching them. 

Interestingly, a clear enhancement during thunderstorms was recently detected by neu- 
tron monitors installed at Mt. Aragats at an altitude of 3250 m a.s.l. in Armenia 
In addition, plastic scintillator-based detectors, arranged close to the neutron monitors, de- 
tected long-duration (10—20 min) 7 rays extending to 40—50 MeV. Generally, owing to 
its detection method, a neutron monitor is believed to be very sensitive to nucleons but 
insensitive to 7 rays and electrons. Thus, Chilingarian et al. [5] have concluded that the 
observed increase of the Aragats neutron monitor is attributable to neutrons generated via 
the photonuclear reaction. 

Similar to the Armenia case, clear enhancements were occasionally obtained by some 
detectors installed at the Yangbajing Cosmic Ray Observatory (30.103°N, 90.523°E; cutoff 
rigidity = 14 GV), which is located on a mountain 4.3 km a.s.l. in Tibet, China. Actually, 
because two electric-field mills (BOLTEK EFM-100) were installed in February 2010 at the 
observatory, five large count enhancements were found to be associated with electric-field 
variations in the rainy season. In this paper we present one prolonged count increase with 
duration of 40 min, which obtained by both the Yangbajing neutron monitor (YBJ NM) and 
a solar neutron telescope (SNT). Utilizing the event, we especially discuss how the observed 
signals are attributed to 7 rays and neutrons produced via the photonuclear reaction. Then, 
we deduce fluxes of 7 rays and neutrons, and compare them with those obtained from other 
experiments and Monte Carlo predictions. 



II. EXPERIMENT 



Due to its high altitude (4300 m a.s.l.) and meteorological conditions from May to 
October, the sky above the Yangbajing Cosmic Ray Observatory is frequently covered with 
thunderclouds. The observatory has three independent detectors: the Tibet air shower 
array 19], YBJ NM 20], and SNT 2l| . The air shower array, working successfully since 
1990, mainly observes 10 12 — 10 17 eV primary cosmic rays. On the other hand, YBJ NM 
and SNT have been operating since 1998, mainly aiming at detecting >100 MeV solar-flare 
neutrons and protons to elucidate the ion-acceleration mechanism in solar flares. YBJ NM 
and SNT are placed close to each other in one building. 
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A. Yangbajing neutron monitor 



YBJ NM consists of 28 NM64-type detectors 22|, |23J having the largest area of 32 m 2 
among world-wide neutron monitors. An NM64 neutron monitor is composed of a BF3 
counter, which is surrounded by polyethylene [(C 2 H 4 ) n ] plates of thickness 7.5 cm and lead 
blocks with an average thickness of 120 gem -2 . The polyethylene plates reflect low-energy 
nucleons accidentally produced in substances close to the detector, while the lead blocks 
multiply impinging neutrons via inelastic scattering processes. 

Each BF 3 counter that contains the BF 3 gas with the density of 3 x 10~ 4 gem -3 has a 
length of 190.8 cm and radius of 7.4 cm. The counter can easily detect a thermal neutron via 
a neutron capture reaction as 10 B + x n — > 4 He + 7 Li, because the cross section of the capture 
reaction increases rapidly as the kinetic energy of the neutron decreases to thermal energy. 
To efficiently decelerate neutrons to thermal energy by elastic collisions with hydrogen nuclei, 
each BF 3 counter is inserted into an additional polyethylene tube with a thickness of 2 cm. 
A 4 He ion created by a neutron capture reaction produces a large amount of ionization loss 
by ~1 MeV or higher in the BF3 counters to provide a sufficiently large signal on its anode. 
Due to the multiplication and thermalization of the incident neutron, the large signal has no 
information about the incident energy. However, the signal can be easily distinguished from 
charged secondary cosmic-ray background events (mainly muons), which provide a small 
signal of ~9 keV. Output signals from individual counters are fed to the data acquisition 
system, and the event number of individual counters is recorded every second. 

It is widely believed that neutron monitors have no sensitivity to electromagnetic com- 
ponents because of the thick lead blocks. However, a photonuclear reaction between 7 rays 
and lead nuclei begins at the 7-ray energy of 7 MeV, and peaking at ~13 MeV [16|. Thus, 
high-energy 7 rays associated with energies >7 MeV can produce neutrons via the photonu- 
clear reaction. Accordingly, neutron monitors might capture such photonuclear neutrons 
produced by thundercloud-related 7 rays extending to 10 MeV or higher energies. 

In order to investigate this possibility, we performed a Monte Carlo simulation based on 



GEANT4 



24j and derived detection efficiencies of an NM64 neutron monitor, including YBJ 



NM, for neutrons, 7 rays, electrons, and positrons in a wide energy range of 1 keV— 1 GeV. 



For this purpose, a geometry of a standard NM64 neutron monitor 23J was constructed, 
and 10 5 mono-energetic particles for each species were illuminated on the same area as the 



neutron monitor. In one mono-energetic simulation, an irradiated particle was isotropically 
injected toward the neutron monitor from the vertical direction to 60 degrees. We choose 
in each simulation (including air propagation simulations described later) a hadronic model 
of QGSP_BERT_HP provided by GEANT4 to treat physics processes of neutrons in the 
atmosphere. 

Figure [TJ shows detection efficiencies determined in this manner for the four particles. 
The present efficiency for neutrons (black circles) agrees well with that obtained by another 
detector simulation conducted by Clem and Dorman 25( (dashed lines). The difference 
in efficiencies at 10 MeV— 1 GeV of the two simulations is maximum 30%. In addition, 
our results for neutrons can reproduce well efficiencies experimentally determined using an 
accelerator neutron beam 26J]. These consistencies validate our simulation results. 

As expected, the present simulation reveals that an NM64 neutron monitor has sensi- 
tivity to electromagnetic components in energies range 10 MeV— 1 GeV. Compared with 
the efficiencies for neutrons, those for 7 rays in the energy range are lower by a factor of 
1/125—1/20. Similarly, high-energy electrons and positrons entering the lead blocks emit 7 
rays via bremsstrahlung, which in turn generate either neutrons via the photonuclear process 
or electrons via pair creation. Since the critical energy of electrons in lead is 7 MeV, these 
cascading processes would continue until energies of electrons and 7 rays are below the crit- 
ical energy. As the incident energy of electromagnetic components increases, the cascading 
becomes more effective in causing the photonuclear reaction. Thus, detection efficiencies for 
electromagnetic components increase (Fig. [JJ. 



B. Tibet solar neutron telescope 

Here, we provide minimal information necessary to understand events reported in this 
paper; detailed information on the Tibet SNT, including detection efficiencies for neutrons 



and 7 rays, is presented in Muraki et al. 



2l| 



SNT installed at the observatory is part of the international solar-neutron observation 
network. It is composed of nine plastic scintillation counters and proportional counters that 
are placed around them. A plastic scintillation counter contains plastic scintillator blocks 
of area and thickness of 1 m 2 and 40 cm, respectively. Thus, the total area of the plastic 
scintillators is 9 m 2 . The counter has a 012.7 cm photomultiplier at the top of the counter 
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for collecting light emissions originating from incident particles. 

Incident charged particles deposit their energies in the thick plastic scintillators via ion- 
ization loss, and hence can be readily observed with SNT. Incident neutrons produce recoil 
ions by scattering protons or carbons in the plastic scintillators, while 7 rays produce elec- 
trons via Compton scattering or pair creation. Through these processes, SNT is able to 
measure neutrons and 7 rays, although it does not differentiate between them. In addition, 
output signals from the photomultiplier are fed to the data acquisition system, amplified and 
discriminated at 4 levels, which correspond to energy deposits of an incident particle of >40, 
>80, >120, and >160 MeV. For each of the nine plastic scintillation counters, individual 
discriminated logical signal is counted by scalers every second. 

Proportional counters complement the plastic scintillation counters. A proportional 
counter has a length of 330 cm and radius of 5 cm, and contains 90% Ar and 10% CH 4 . 
Thirty proportional counters are placed above the nine plastic scintillation counters, while 
seventy-two proportional ones shield the 4 sides of the plastic scintillation ones. Therefore, 
the surrounding counters can be utilized as an anti-counter to separate photons and neu- 
trons from charged particles. In fact, using the surrounding proportional counter signals 
in anti-coincidnece, the four discriminated counting rates of the central plastic counters are 
reduced by a factor of 0.2—0.25. 

C. Electric-field mill 

To measure electric-field variations, two commercial electric-field mills (BOLTEK EFM- 
100) were installed on the premises. One is mounted on the ground, while the other is located 
on the roof of a central building; hereafter, denoted as EFM1 and EFM2, respectively. The 
two mills are arranged ~25 m apart with a vertical distance of 3.4 m. Individual output 
signals are transmitted to the central building with optical cables, directly fed to PCs, and 
recorded every 0.1 s as the electric field strength in the range ±40 kVm -1 with a resolution 
of 20 Vm" 1 . 

The electric field strength measured by EFM2 is always higher, by a factor of ~2, than 
that by EFM1. Such an enhancement of an electric field is often caused by distortion of 
local electric field lines because of obstructions such as a building. In fact, EFM2 is installed 
near a corner of the roof of a building, and hence are more largely affected by such distortion 
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than EFM1, which is located on the ground with few surrounding obstructions. Considering 
this disparity, we use only EFM1 data in this paper. 



III. OBSERVATIONS 



overview 



Examining the data over the 2010 rainy season from May to October, we visually found 
25 events in which electric fields largely deviate from fair-weather states of <100 Vm" 1 . 
Five of them accompanied prolonged count enhancements. Three of the five events are 
clearly observed by either YBJ NM or SNT, lasting for 10—20 minutes. Similar events with 
duration of 10—20 minutes have been already reported by other measurements (for example, 
3|, |5J). On the other hand, the remaining two events, detected by both YBJ NM and SNT, 
last for >30 minutes. Such a long-lasting emission has never been observed. 

With two reasons, we selected one of the two events that are detected by both YBJ 
NM and SNT. One is that the event is a fast observation of the longest-duration emission 
among other long-duration events. The other is that the selected event clearly correlate 
with electric fields (as shown later), while the other has only a poor correlation with electric 
fields. Although statistical significance of the latter event for YBJ NM and SNT(>40 MeV) 
were around two times higher than the selected event, we will have to collect additional ones 
in order to well understand the nature of such a poorly correlated event. 

B. Count histories 

Figure [2] shows five-minute counting rates by YBJ NM and SNT and one-second electric- 
field variations obtained over 3:00—7:00 UT on July 22, 2010. All counting rates by YBJ 
NM and SNT are corrected for atmospheric pressure variations. In this event, YBJ NM 
count rates and SNT ones in >40 MeV clearly increase around 5 UT [Figs. [2](a) and (b)]. 
In addition, higher-energy channels of SNT [panels (b) — (e) of Figs. [2] appear to show count 
enhancements in coincidence with the above clear increases. Given the clear signals, in 
particular, the >40 MeV channel of SNT vetoing charged particles with the anti-counter 
[Fig. [2^b)], we can conclude that >40 MeV 7 rays and/or neutrons reach the detectors to 
produce the observed signals. 
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With a criterion that individual counts of the >40 MeV channel of SNT continuously 
have >2<j statistical significance above background, we define burst time as 40 min at 
4:30—5:10 UT. Here, by excluding the data in this period and fitting the remaining data 
with a quadratic function, we estimate the background for YBJ NM and SNT (gray dashed 
curves in Figs. [2]). Subtracting the interpolated background from total observed counts in 
the burst period, we obtain net count increases for the burst recorded by YBJ NM and SNT; 
these are listed in Table [I] together with their statistical significance. Hereafter, the burst is 
simply called 100722. 

Generally, the counting rate of a neutron monitor, including YBJ NM, does not sim- 
ply obey Poisson distribution because of the multiplication of one incident neutron in the 
lead blocks. Usoskin et al. {2?! provide a more detailed explanation on how these effects 
cause non-Poissonin fluctuations in NM data. They described that a statistical significance 
obtained by a NM usually should be reduced by a factor of 1.2—2, depending on the geo- 
magnetic cut-off rigidity and the atmospheric depth at NM locations. Thus, the statistical 
significance obtained (Table [I]) may decrease by half. Importantly, both YBJ NM and SNT 
simultaneously recorded large enhancements in association with electric-field variations. 

Based on the following features of the event observed, we may conclude that it is asso- 
ciated with thunderclouds, but not lightning. First is its long duration; apparently, such 
a long-duration emission would not be generated by lightning and/or its related phenom- 
ena that generally last for milliseconds or less. Second is that the electric field strength in 
the burst period does not change rapidly (within 1 s), but gradually [Figs. EJ^e)]. In addi- 



tion, a 
7 



though not homogeneous, these features have already been reported by many groups 



11 



as thundercloud-related emissions. 



C. Relation with electric fields 

Figure [3] represents detailed time variations of YBJ NM and SNT, together with the 
averaged electric fields. Clearly, peaks of YBJ NM and SNT signals for 100722, obtained 
over 4:50—5:04 UT, correlate with those of electric fields in the same interval [Figs. EJ^a) and 
(b)]. For clarity, Figure H] shows the correlation between the present burst and the electric 
field measured by EFM1. We computed a correlation coefficient between the count variations 
of YBJ NM and SNT and the electric field as 0.79 (0.01) and 0.77 (0.03), respectively. Each 
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number in parentheses represents a correlation coefficient outside the burst period. 

An electric field in the downward direction is measured as a positive field. Thus, the 
positive electric fields correspond to the existence of positive charges overhead, which are 

when thunderclouds exist at a mature 



frequently observed at Tibet [33|] and New Mexico 
stage over a field mill on the ground. Furthermore, such a thundercloud generally forms 
tripole electrical structures, which consist of positive, negative, and positive layers from top 
to bottom, which in turn accelerates electrons therein toward the ground. 



IV. NEUTRON PRODUCTION AND PROPAGATION IN THE AIR 
A. outline 

According to Babich et al. {35 1, a yield rate of a photonuclear neutron per one gamma 
ray with energies >10 MeV is 4.3 x 10~ 3 . Produced neutrons propagate in the atmosphere, 
attenuated by elastic and/or inelastic collisions with air nulcei. Assuming neutrons propa- 
gate over L =1 km (0.1 km) to reach the observatory, neutrons produced decrease in number 
by a factor of exp (— L/X Q ) = 2 x 10~ 3 (0.5). Here, A n represents an attenuation length of 
neutrons in the atmosphere, calculated as A n = 13 gem -2 for 20-MeV neutrons using total 
cross section between a neutron and an air nucleus 36J. As a result, a >10-MeV 7 ray is 
found to produce 10~ 5 — 10 -3 neutrons to arrive at the observatory. Given this arrival rate of 
neutrons and derived detections efficiencies for neutrons and 7 rays (Fig. [1]), we expect that 
>10-MeV 7 rays would be able to considerably contribute to the signals detected by YBJ 
NM. To better understand how much photonuclear neutrons propagate to the observatory, 
we performed a GEANT4 simulation. 

For the purpose of simulating neutron production via the photonuclear reaction and 
neutron propagation in the atmosphere, we constructed five atmospheric layers starting 
from the observatory level (4.3 km a.s.l.) to 5 km higher. Each rectangular atmospheric 
layer has a vertical length (z direction) of 1 km and horizontal length (xy directions) of 10 
km, and consists of N 2 , 2 , and Ar with mole ratios of 78.1%, 21.0%, and 0.9%, respectively. 
Air density in the individual layers is fixed at 7.7 x 10~ 4 gem -3 for 4.3—5.3 km a.s.l., 7.0 x 
10 -4 gem -3 for 5.3—6.3 km a.s.l., 6.6 x 10 -4 gem -3 for 6.3—7.3 km a.s.l., 5.6 x 10 -4 gem -3 
for 7.3—8.3 km a.s.l., and 5.0 x 10 -4 gem -3 for 8.3—9.3 km a.s.l [37]]. In the following 
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simulations, seven source heights are assumed to be 0.3, 0.6, 0.9, 1.5, 2, 3, and 5 km 
above the observatory level. From each source height, one million 7 rays were injected to 
the atmosphere to produce secondary particles. The secondary particles, propagating to the 
observatory, were saved with their species, energy, x-y positions, azimuth, and zenith angles. 



Bremsstrahlung 7 rays derived from runaway electrons have been thoug 
exponentially cut-off power-law spectrum, with a cut-off energy of ~7 MeV 



it to have an 
. How- 



38 



ever, the recent AGILE observation 



31] indicated that a high-energy part (>1 MeV) of 



the TGF spectrum extending from 10 MeV to 100 MeV can be explained by a power-law 
spectrum with a spectrum index, 0, of —2.7 rather than an exponentially cut-off one. Sea- 
level observations of long-duration 7 rays also showed that a source 7-ray spectrum may be 
described as a power-law type with ~ —2 Theoretically, of a bremsstrahlung 7-ray 
spectrum has the hardest limit of —1. We therefore assumed a power- law spectrum as an 
initial photon spectrum in this study and is — 1, —2, or —3. The minimum and maximum 
energies of the spectrum are set at 10 and 300 MeV, respectively, to fully cover the presently 
relevant energy range. In addition, downward directions of initial 7 rays were assumed to 
be distributed either isotropically within 0—30 degrees or over a Gaussian beam with a half- 
opening angle of 30 deg. Both types would be expected from runaway electrons moving in 
electric fields in air, because moving electrons are subjected to multiple scatterings with air 
molecules, and the geo metrical or electrical structure of electric fields in thunderclouds may 
not be very simple 40]. 



B. Energy spectrum 

Figures and E] show neutron energy spectra obtained by the isotropic and Gaussian 
angular distributions, respectively. There is no significant difference in shape of the neutron 
spectra between the two angular distributions. These neutron spectra suggest that the 
neutrons arriving at the observatory have a mean energy of 1—10 MeV and the maximum 
energy of produced neutrons is about one-third of that of the 7 rays emitted from a source. 
The former feature has been reported by Carlson et al. and Babich et al. 35] as well. 

Figure H represents spectra for 7-ray and electrons assuming the isotropic emission of 
initial 7 rays. Similar to neutron spectra, those for 7 rays and electrons, assuming the 
Gaussian beam emission, do not largely change from the isotropic ones. 
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C. Survival rate 



Figure [H] shows survival rates for >1 keV neutrons and >10 MeV 7 rays for the two 
angular distributions, sampling arriving neutrons (7 rays) with energies of >1 keV (>10 
MeV) and calculating a ratio of the number of the arriving neutrons (7 rays) to that of 
primary 7 rays. The threshold energy of 1 keV for neutrons does not affect our results, 
because neutrons with energies <1 keV constitute a maximum 5% of all neutrons produced. 
As expected, the neutron survival rates for the two angular distributions are similar in shape 
and intensity, having at most 10% difference in rate. Depending on spectrum indices, the 
neutron survival rates are generally constant at ~10 -3 until the source height is around 1 
km, and then decrease to 10~ 4 — 10 -5 . The derived survival rates quite agree with those 
simply calculated in Sec. IIV Al 

As can be easily seen, each neutron survival rate has its peak at the source height of ~0.6 
km, which corresponds to ~50 gem" 2 . The shape of the survival rates of neutrons simply 
reflects the product of the probability that the photonuclear reaction occurs at the point 7 
rays propagate in the atmosphere and that the produced neutrons are attenuated, which is 
proportional to [1 — exp (— H/\)] x exp (—H/X n ). Here, H represents the assumed source 
height, while A p and A n represent the interaction length of 7 rays to cause photonuclear 
reaction, which is ~3000 gcm~ 2 at the peak cross section of 15 mb, and the attenuation 
length of neutrons, respectively. A n in the relevant neutron energies of 1—100 MeV range 
between 20—100 gem -2 , corresponding to 0.2—1.4 km. 

V. CONTRIBUTION RATIOS TO THE SIGNALS 
A. Method 

Given the simulated neutron spectra (Figs. |5]and|6]), and those of 7 rays and electrons 
(Fig. [7]), as well as the detection efficiencies of the neutron monitor (Fig. [TJ), we can examine 
how neutrons and electromagnetic components contribute to signals that are expected to be 
detected by YBJ NM and SNT during thunderstorms. 

As argued so far, we presume that the four components, neutrons, 7 rays, electrons, and 
positrons, explain the count increases observed by YBJ NM, and that neutrons and 7 rays 
contribute to SNT signals because SNT utilizes the anti-counter to reject charged particles. 
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rK 2 

Am(t) = aA / I^SiiE^E^dEi, (1) 



Therefore, a predicted count increase at a given time t for individual particles, Arii(t), is 
written by 

«if a 

assuming that the relevant particles have the same production history, I(t), to be generated 
during thunderstorms. Here, a represents a normalization factor with the unit MeV _1 s _1 m -2 
for a source spectrum and A represents the area of YBJ NM (32 m 2 ) or SNT (9 m 2 ). Ei 
denotes the energy of a particle type i, Si(Ei) represents the spectra (Fig. [5] for isotropic 
emissions), and €{(Ei) denotes the detection efficiencies of YBJ NM (Fig. [T]) or SNT (Fig. 5 



of |2l|). In the present study, K% is set to 1 keV for neutrons and 10 MeV for electromagnetic 
components, while K 2 is fixed at 300 MeV for all components. By integrating Eq.flTJ) over 
a certain time interval of t2—ti, we can obtain an expected net count increase due to each 



particle (ANA as 



AN = I Am(t)dt. (2) 



Under the present assumption, the simulated spectra Si(E{) in Eq. (JTJ are independent 
of time t. Accordingly, a ratio of AN/ AN is calculated as 



AN Iki Si(Ei)ei(Ei)dEi 



SkI Si(Ei)ei(Ei)dEi ' 
which shows the contribution fraction of each species to an expected signal. 



B. YBJ NM signals 

Figure [9] depicts contribution ratios of neutrons and 7 rays for YBJ NM, assuming an 
isotropic angular distribution. As expected, contribution ratios for the Gaussian angular 
distribution are almost the same as those for the isotropic distribution. For clarity, con- 
tribution ratios for electrons and positrons are not shown. Interestingly, the contribution 
ratios of neutrons and 7 rays do not depend largely on (3. Therefore, it is obvious that 7 
rays dominate (96% to 85%) the fraction of the expected count increase as the source is 
farther, while neutrons contribute a maximum of 15%. 
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C. SNT signals 



Similarly, contribution ratios for SNT signals can be calculated using detection efficiencies 
for neutrons and 7 rays in Fig. 5 of Muraki et al. [2l| . The Ninety-nine percent of the 
observed signal for >40 MeV channel of SNT is dominated by 7 rays, while the remaining 
three higher energy channels are almost fully contributed by 7 rays. These results for SNT 
are mainly ascribed to a relatively small fraction (<5%) of neutrons produced in >40 MeV 
energies via the photonuclear reaction (Figs. and E]) • 

VI. TIME HISTORIES OF YBJ NM AND SNT 

I(t) can be naturally assumed to follow the time history of the electric field [Fig. EJ^f)]. 
Utilizing the one-sec electric-field variations as I{t), Eq.([T]) has only one unknown parameter, 
a, that needs to be determined by comparing an expected time profile of YBJ NM or SNT 
to the observed one- minute profile. Further, we test the following two hypotheses. First, 
the relevant particles are produced only when the electric field at the surface has positive 
polarity and electrons in thunderclouds are accelerated toward the ground. Second, the 
particles are generated when the field has negative and positive polarities. 

For the purpose of introducing the mathematical form of I(t) for the first assumption, 
the positive electric field strength [Fig. I21^f)] in the burst periods is divided by the maxi- 
mum strength of 26.8 kVm -1 , and the negative electric field strength and that outside the 
individual burst periods are set to zero. On the other hand, absolute values of the electric 
field strength divided by the above mentioned maxima, are considered as I(t) for the second 
assumption. In addition, I(t) in this case is zero outside the burst times. Therefore, J I(t)dt 
for both assumptions is normalized to one. Hereafter, we call the first and the second as- 
sumptions "negative emission" and "bipolar emission", respectively. Substituting each I(t) 
in Eq. ([TJ and integrating Eq. (P2J) every 60 s in individual burst intervals, we can prepare a 
one- minute expected time profile depending on each I(t) and compare it with the observed 
one-minute counting rate of YBJ NM and SNT (Fig. [3]). Next, we compute 



to search for the \ 2 minimum with a being a free parameter. Here, iV b(ti) and N ex (ti) 
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represent the observed one-minute counts of YBJ NM or SNT and the model-predicted 
counts at a given time ti in the burst intervals, respectively. Statistical errors associated 
with iV ob (£j) are written as 0" o b(tj). Summation was carried out over each burst interval. 

In fact, each I(t) produces the same shape of predicted count history and the same \ 2 
minimum for YBJ NM or SNT, despite using simulated spectra [S{ in Eq.flTJ] obtained with 
various sets of f3 and H. This is because Eq. ([I]) has only one unknown parameter a, and the 
shape of Si in Eq. (JTJ is independent of t. To concretely determine a, we first independently 
evaluated a for YBJ NM and SNT with each I(t) by the above method, using simulated 
spectra obtained by 21 combinations of ((3, H); (3 of —1, —2, —3 and H of 0.3, 0.6, 0.9, 
1.5, 2, 3, and 5 km. The derived a for YBJ NM and SNT are shown in Figure [TU1 Next, 
subtracting a acquired from YBJ NM data with a set of (/3, H) from that acquired from 
SNT data with the same set of (/3, H), we searched for the smallest difference in a obtained 
by the two independent detectors. 

As clearly seen in Figure [101 a difference in a for YBJ NM and SNT is the smallest 
at (f3, H) = (—2, 900 m) and (—2, 600 m) for the negative emission and bipolar emission, 
respectively. Table [III displays the calculated a and \ 2 minima, together with a set of /3 
and H. Figures [11] compare the observed count histories of YBJ NM and SNT with those 
expected from the parameters listed in Table [TT] under the two emissions. The x 2 values 
(Table [IT]) clearly suggest that the observed time profiles are reproduced by the negative 
emission rather than the bipolar emission. 



VII. DISCUSSION 



A. Gamma-ray emissions 

1. Characteristics of 7 -ray emissions 

The present study revealed that high-energy 7 rays with energies >40 MeV originate from 
summer thunderstorms. According to sea- level observations in winter thunderstorms 
long-duration 7 ray emissions from winter thunderstorms extend to 10—20 MeV. This may 
be due to a difference in atmospheric density at grou nd and high mountains. In fact, a TGF 
spectrum averagely extends to few tens of MeV [2^, [jo]], or 100 MeV on rare occasions {31) ]. 
It is believed that TGFs occur at altitudes of 15—20 km {32]. These results including the 
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present one may imply that a lower atmospheric density is attributable to a higher 7-ray 
emission. 

Compared with other thundercloud-related 7-ray events, the duration of the present 
event is exceptional with its long duration of 40 min. Wash-out radioactive radon and its 
decay products frequently cause count increases in ground-based detectors. In addition, a 
duration of such a radon effect is around 20—30 min corresponding to their half-lives. Thus, 
the duration of the radon effect is similar to the present one. However, the radon families 
generate <3 MeV 7 rays, being unable to give signals in YBJ NM and SNT. 



According to electric field measurement at the Tibet plateau (4.5 km a.s.l.) 33] and a 
mountain in New Mexico (3.2 km a.s.l.) [34], the mature stage of summer thunderclouds 
seems to last for ~1 h. In addition, the measurement in New Mexico revealed that a vertical 
potential relative to the surface in the mature stage is quasi stable which are required for 
electrons to be continuously accelerated in thunderclouds in order to produce prolonged 
7-ray emissions. Therefore, we infer that the present event is mainly associated with the 
mature stage of the Yangbajing thunderstorms. On the other hand, mature stages of winter 
thunderstorms at a costal area of Japan sea last for <10 min [4l|]. In fact^all thundercloud- 
related 7 rays observed in winter lasted for at most a few minutes 0,3, S]- Thus, it is 
deduced that the longevity of the mature stage plays an important role in determining the 
duration of thundercloud-related 7 rays. 

From the 7-ray emission of the present event, a source height H was estimated as H = 900 
m (Tabled]), giving the source altitude of 5.2 km a.s.l. Qie et al. {33]] reported that a cloud 
base of summer thunderclouds above the Tibetan plateau (4.5 km a.s.l.) is generally located 



at ~1 km. In addition, Marshall et al. [42| clearly showed that a bottom positive layer of 
a summer thundercloud in New Mexico is located at 4.5—5.5 km a.s.l. Thus, the source 
altitude of 5.2 km a.s.l. is in good agreement with altitudes of the cloud base and the 
positive bottom layer obtained from these observations. 

As clearly seen in Figure [2l time structures for YBJ NM and SNT are different with 
each other. In particular, YBJ NM showed no count increases at the burst onset, while 
all the SNT channels (>40 MeV to >160 MeV) provided count enhancements in 5—10 
minutes after the onset. These peculiar time structures might be caused by moving of 
thunderclouds and limited illumination of higher-energy part of bremsstrahlung gamma rays 
emitted from thunderclouds. Actually, it is confirmed that long-duration gamma rays move 
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with thunderclouds |7j. In addition, the bremsstrahlung gamma rays, especially gamma rays 
with an energy being close to that of accelerated electrons, would be relativistically beamed 
into a narrow cone. For example, a half opening angle of the cone, 9 ~ I/7, is 0.1° for a 300 
MeV electron, where 7 is Lorentz factor. Given H = 900 m, we can obtain a radius of the 
gamma rays arriving at the observatory as at most 1.6 m (900 mx tan 0.1°). Because YBJ 
NM is located ~10 m apart from SNT, >40 MeV gamma rays moving with thunderclouds 
might not happen to face towrad YBJ NM over the burst onset. 



2. Electric potential 

Due to ionization loss of electrons, an electric potential of 40 MV is not high enough to 
accelerate electrons to 40 MeV. In practice, an electric field strength of 240—270 kVm^ 1 is 
required for electrons of 1 — 10 MeV to be accelerated to 40 MeV assuming a vertical length 



42 



43j. 



of a high-electric field region is 0.5 — 1 km, as determined by balloon experiments 
Multiplying this field strength by the assumed vertical length, the electric potential of at least 
120 MV must be established in the thunderclouds. This value of 120 MV is approximately 



equal to the maximum potential of 130 MV observed by balloon soundings 44j. In addition, 
the AGILE observation of TGFs showed that the electric potential in thunderstorms is on 
the order of 100 MV [3jJ over macroscopic lengths such as cloud sizes or intracloud distances. 
Accordingly, the present observations may show manifestation of the highest potential field 
during thunderstorms. 



3. Avalanche multiplication factor 

In addition to quasi-stable electric fields, a stable or quasi-stable source of seed electrons 



would generally be needed for prolonged 7-ray emissions. Gurevich et al. js] originally 
postulated that secondary cosmic rays consist of seed electrons, which increase in number 
and emit bremsstrahlung 7 rays. Thus, according to this premise, we derive an avalanche 
multiplication factor, M, expected from the RREA mechanism. 

Using a and (3 for the negative emission (Table HI!), a source 7-ray spectrum, F S (E), 
can be described as F S (E) = a w E~ 2 . Here, E is a photon energy in MeV and a w = 
(4.3 ± 0.2) x 10 3 m _2 s _1 MeV _1 is a weighted mean calculated by the values of a from YBJ 
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NM and SNT. Using F s and the burst duration AT = 2400 s, we estimated the total number 
of electrons with 10—100 MeV energies as N e ~ 10 14 , in the same manner as estimated 
by Tsuchiya et al. [7|. For this purpose we assumed a single acceleration region in the 
thundercloud with the vertical length and horizontal one of Z = 500 m or 1000 m and L = 
600 m [7], respectively. In reality, a positive or a negative charge layer of thunderclouds 
may consist of multi cells [e.g. 4^] to form several acceleration area therein. Thus, the single 
acceleration region is a simple assumption to consider individual particle accelerations. 

The secondary cosmic-ray electron flux above 1 MeV at the relevant altitude is Iq ~ 400 
m -2 s _1 |46|. Therefore, the number of such electrons N cr entering the acceleration region 
in the burst period is computed as 



N cr = I x A s x AT ~3x 10 
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L 



600 m 



giving M as 



M = N C /N CI = 300 



600 m N 2 



L 

Furthermore, based on the RREA mechanism, M thus derived is described as 

/r7m , 7300 kV 

M = exp(ZX),X= — -, 3 

PV 1 h V- 276 P/l atm V ; 



where A and V denote a length parameter given by Dwyer [10[ and electric field strength 

in kVm" 1 , respectively. Substituting M = 300, Z = 500 or 1000 m, and P = 0.55 atm 

(average pressure at H = 900 m) in Eq. 03]), we obtain V = 240 and 190 kVm -1 for Z = 500 

m and 1000 m, respectively. These values of V are consistent with the above estimated field 

strength to accelerate electrons to 40 MeV or higher energies. 

Conducting sea-level observations in winter, Tsuchiya et al. |7j showed that secondary 

cosmic-ray electrons are multiplied by a factor of 3—30 to produce thundercloud-related 7 

I 

rays. On the other hand, Chilingarian et al. [5[ obtained a multiplication factor of ~330 
with a high- mountain measurement in summer. From these results as well as our result, a 
multiplication factor in high mountains can be considered to be different from that at sea 
level. However, the above M becomes 30 if L = 2 km, which is observed as the horizontal 



extent of a bottom positive layer in a summer thundercloud 42] . Thus, if L is longer than 2 



km, the estimated M may become consistent with that derived from sea-level observations. 
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B. Neutron emissions 



1. Comparison with the Aragats neutron monitor 

Similar to the present event, Chilingarian et al. jlj] demonstrated that a count enhance- 
ment lasting for 10 min was detected by the Aragats neutron monitor located at 3250 m 
a.s.l. As a result, they concluded that the observed increase is fully attributable to neutrons 
related to the photonuclear reaction. On the other hand, the present results demonstrate 
that >10-MeV 7 rays dominate the signals observed by YBJ NM. This is a main difference 
between the present study and that by Chilingarian et al. |5j. 

The present simulation clearly showed that an NM64 neutron monitor, which was also 
used by Chilingarian et al. |5], has low, but not negligible, sensitivity to 7 rays. In addition, 
the survival probability of neutrons and 7 rays at the Aragats observatory would not largely 
change from the present one (Fig. [8]), because the air density at the Aragats observatory, 
which is ~9 x 10 -4 gem -3 , is not very different from that at the Yangbajing site, which 
is ~8 x 10~ 4 gem -3 . In fact, using the GEANT4 simulation, Chilingarian et al. [5| de- 
rived 2.3 x 10 -3 as survival probability of neutrons arriving at their observatory, assuming 
bremsstrahlung 7 rays propagate over 1500 m. This value is nearly consistent with the sur- 
vival probability of neutrons that is derived in the present study for H = 1500 m, which is 
5 x 10~ 4 — 2 x 10 -3 (Fig. [8]). Consequently, not neutrons but 7 rays may possibly dominate 
enhancements detected by the Aragats neutron monitor. 



2. Number of neutrons produced 

Using the derived value of a w , we evaluate the fluence of neutrons, f n , arriving at the 
observatory in energies 1 keV— 300 MeV, by the following formula: 

p /-300 MeV 

f n = a w AT / I(t)dt / S n (E Q )dE Q = 1.4 x 10 4 m" 2 , 

J AT Jl keV 

where AT = 2400 s and S n (E Q ) represents the simulated neutron spectrum, assuming (3 
and H are —2 and 900 m (Table [II]), respectively. Carlson et al. 17] and Babich et al. 
H ™ ft* p—e ar _ produced by energetic 7 rays are observable at 
ground level when a 7-ray source is locate <5 km, since the neutron fluence is expected as 
(0.03— l)xl0 4 m -2 for the former prediction and 10 3 — 10 7 m~ 2 for the latter one. Actually, 
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the value of f n is consistent with their predictions. Thus, this agreement imply that the 
photonuclear reaction certainly occurs during mature stages of thunderclouds. 



VIII. SUMMARY 



The prolonged 7-ray event, lasting for 40 min, was observed on 2010 July 22 at Yangbajing 
in Tibet, China. Such a long-duration event associated with thunderstorms have never been 
observed. In addition, the present observations clearly showed that 7 ray extending to 
energies >40 MeV were detected by SNT and very likely by YBJ NM. Given these results, 
the present emissions strongly suggest that electrons are accelerated beyond at least 40 
MeV in 40 min, by quasi-stable electric fields, which were formed during the mature stage of 
summer thunderclouds. The present duration is at least 5 times longer than those observed 
in winter thunderstorms at the coastal area of the Japan sea. Probably, one of the main 
reasons for this difference would be ascribed to a difference in life cycles of mature stages of 
winter and summer thunderclouds. 

The high-energy 7 rays would produce neutrons via the photonuclear reaction of 14 N(7, n) 13 N. 
The present simulation showed that the arriving neutron flux at >1 keV is expected to be 
lower than that of arriving 7 rays at >10 MeV by more than two orders of magnitude. 
Moreover, it revealed that unlike previously believed, neutron monitors are not insensitive 
to 7 rays. Consequently, it is found that bremsstrahlung 7 rays largely attribute the signal 
obtained by YBJ NM and photonuclear neutrons give only a small contribution to the 



signal. The present study demonstrated that world-wide networks of neutron monitors 
and solar neutron telescopes 
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are useful for observations of thunderstorm-related 



7-ray emissions. 
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TABLE I. Net count increases and statistical significance. 



A r\Q (significance) 



YBJ NM 34000 ± 4200 (8.1a) 

SNT >40 MeV 44000 ± 3500 (13a) 

SNT >80 MeV 16000 ± 2400 (6.7a) 

SNT >120 MeV 8700 ± 1500 (5.8a) 

SNT >160 MeV 4600 ± 970 (4.7a) 

a Each quoted error includes fluctuations of the background and total observed counts. 
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TABLE II. x 2 minima and spectrum parameters determined. 



Negative emission Bipolar emission 



YBJ NM SNT YBJ NM SNT 

X 2 /d.o.f. 49.4/39 46.2/39 110/39 59.7/39 

a (xlO 3 MeV _1 m- 2 s -1 ^ 4.3 ± 0.3 4.2 ± 0.3 2.7 ±0.2 2.8 ± 0.2 

(/£,H1) (-2,0.9) (-2,0.6) 



a A normalization factor of an assumed power-law gamma-ray spectrum. 
b An estimated photon index of a power-law gamma-ray spectrum. 
c A source height (km) estimated. 
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FIG. 1. Detection efficiencies of an NM64 neutron monitor for neutrons, 7 rays, electrons, and 
positrons, as determined by the GEANT4 simulation. A dashed curve indicates detection efficiency 



of an NM64 neutron monitor obtained by Clem and Dorman 25f] . Open circles represent experi 



mental results by Shibata et al. 
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FIG. 2. Count rates per 5 min detected by YBJ NM and SNT, and one-second variations of 
EFM1 over 3:00-7:00 UT on July 22, 2010. Panel (a) shows count rates by YBJ NM; while panels 
(b)-(e) show count rates by >40 MeV, >80 MeV, >120 MeV, and >160 MeV SNT with anti- 
coincidence. Panel (f) shows the one-second variations by EFM1. Dashed gray curves in panels 
(a) — (e) indicate the estimated background, while vertical dashed lines in all panels represent the 
defined burst periods. The horizontal axes show universal time. Error bars are statistical la except 
for panel (f). 
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FIG. 3. One-minute count histories in percent observed by YBJ NM and SNT, and one-minute 
averaged electric-field variations. Data points with la error bars in panels (a) and (b) correspond 
to the variations in YBJ NM and >40 MeV SNT with anti-coincidence for 100722, respectively. In 
all panels, histograms (solid lines) represent the average field variations by EFM1. The horizontal 
axes show universal time. Left and right vertical axes denote the count variations in percent and 
the electric field strength in kVm -1 , respectively. Vertical lines in each panel indicate the burst 
onset and end times. 
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FIG. 5. Derived energy spectra of neutrons reaching the observatory level, assuming the isotropic 
emission of initial 7 rays. Left panels show those spectra with an initial source height fixed at a 
given value in each panel. Solid, dashed, and dotted lines in the left panels correspond to (3 of — 1, 
—2, and —3, respectively. Right panels indicate those spectra in which (3 is fixed at a constant 
value with source heights of 0.3 (solid), 1.5 (dashed), and 3 km (dotted). Each horizontal axis 
denotes kinetic energy of neutrons in MeV at the observatory. The vertical axes indicate relative 
values divided by the number of incident 7 rays 1 x 10 6 . 
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FIG. 6. Same as Fig. \5\ but for Gaussian-type angular distribution of initial 7 rays. 
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FIG. 7. Same as Fig. but for gamma rays (left) and electrons (right). 
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FIG. 8. Survival rates of >l-keV neutrons (filled) and >10-MeV 7 rays (open) at the observatory- 
level, obtained by the 7-ray power law spectrum simulations. Left and right panels indicate the 
isotropic and Gaussian beam distributions of initial 7 rays, respectively. Circles, squares, and 
triangles correspond to j3 of —1, —2, and —3, respectively. The horizontal axes show source 
heights. Errors are statistical la errors. 
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FIG. 9. Contribution ratios of 7 rays (squares) and neutrons (circles) for YBJ NM signals, plotted 
against assumed source heights. Left, middle, and right panels correspond to (3 of —1, —2, and —3, 
respectively. Errors are statistical la errors. 
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FIG. 10. Comparison of derived a for YBJ NM (open symbols) and SNT (filled symbols). Left 
panel shows the negative emission, while right one denotes the bipolar emission. Circles, squares, 
and triangles correspond to f3 of —1, —2, and —3, respectively. The horizontal axis in each panel 
represents assumed source height in km. 
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FIG. 11. Comparison of background-subtracted count histories with expected counts for 100722. 
Left and right panels correspond to YBJ NM and SNT (>40 MeV) with anti-coincidence, respec- 
tively. Top and bottom panels represent the negative and bipolar emissions, respectively. Blue 
and red lines in each panel indicate a count history predicted by 7 rays and neutrons, respectively. 
Each magenta line shows an expected count history, which is summed over the counts from the 
relevant particles. Two vertical dashed lines in each panel denote the start and end times of the 
count increase. 
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